# -*- coding: utf-8 -*-+
from PIL import Image
import matplotlib.pyplot as plt
from skimage import transform,io
from scipy import stats,signal,fftpack,ndimage,misc
import numpy as np
import numpy.fft as fp
from skimage import io,color,restoration,img_as_float
import timeit
import matplotlib.pyplot as plt


#1.numpy.fft的频域高斯模糊滤波器
#plt.figure(figsize=(10,11))
#plt.gray()
#im=np.mean(io.imread('F:\\python练习\\图像实战练习\\images\\lena.jpg'),axis=2)
#gauss_kernel=np.outer(signal.gaussian(im.shape[0],5),signal.gaussian(im.shape[1],5))
#freq=fp.fft2(im)
#assert(freq.shape==gauss_kernel.shape)
#freq_kernel=fp.fft2(fp.ifftshift(gauss_kernel))
#convolved=freq*freq_kernel
#im1=fp.ifft2(convolved).real
#plt.subplot(2,3,1),plt.imshow(im),plt.title('Original Image',size=20),plt.axis('off')
#plt.subplot(2,3,2),plt.imshow(gauss_kernel),plt.title('gauss_kernel',size=20),plt.axis('off')
#plt.subplot(2,3,3),plt.imshow(im1),plt.title('Output Image',size=20),plt.axis('off')
#plt.subplot(2,3,4),plt.imshow((20*np.log10(0.1+fp.fftshift(freq))).astype(int))
#plt.title('Orignal Image Spectrum',size=15),plt.axis('off')
#plt.subplot(2,3,5),plt.imshow((20*np.log10(0.1+fp.fftshift(freq_kernel))).astype(int))
#plt.title('gauss_kernel Image',size=15),plt.axis('off')
#plt.subplot(2,3,6),plt.imshow((20*np.log10(0.1+fp.fftshift(convolved))).astype(int))
#plt.title('Output Image Spectrum',size=15),plt.axis('off')
#plt.subplots_adjust(wspace=0.2,hspace=0)
#plt.show()


#二维高斯LPF核频谱
#im=color.rgb2gray(io.imread('F:\\python练习\\图像实战练习\\images\\lena.jpg'))
##高斯
#guass_kernel=np.outer(signal.gaussian(im.shape[0],1),signal.gaussian(im.shape[1],1))
#freq=fftpack.fft2(im)
#freq_kernel=fftpack.fft2(fftpack.ifftshift(guass_kernel))
#plt.imshow((20*np.log10(0.01+fp.fftshift(freq_kernel))).real.astype(int),cmap='coolwarm')
#plt.colorbar()
#plt.show()


#高斯模糊滤波器------signal.fftconvolve()函数
#img=np.mean(io.imread('F:\\python练习\\图像实战练习\\images\\mandrill.jpg'),axis=2)
#guass=np.outer(signal.gaussian(11,3),signal.gaussian(11,3))
# #模糊
#im_blurred=signal.fftconvolve(img,guass,mode='same')
#fig,(ax_orignal,ax_kernel,ax_blurred)=plt.subplots(1,3,figsize=(15,8))
#ax_orignal.imshow(img,cmap='gray')
#ax_orignal.set_title('Original',size=20)
#ax_orignal.set_axis_off()
#ax_kernel.imshow(guass,cmap='gray')
#ax_kernel.set_title('Guassian',size=15)
#ax_kernel.set_axis_off()
#ax_blurred.imshow(im_blurred,cmap='gray')
#ax_blurred.set_title('Blurred',size=20)
#ax_blurred.set_axis_off()


#卷积后绘制原始图像和频谱
#F1=fftpack.fft2((img).astype(float))
#F2=fftpack.fftshift(F1)
#plt.figure(figsize=(15,8))
#plt.subplot(221),plt.imshow((20*np.log10(0.1+F2)).astype(int),cmap=plt.cm.gray)
#plt.title('Origure',size=20)
#plt.axis('off')
#
#F1=fftpack.fft2((im_blurred).astype(float))
#F2=fftpack.fftshift(F1)
#plt.figure(figsize=(15,8))
#plt.subplot(222),plt.imshow((20*np.log10(0.1+F2)).astype(int),cmap=plt.cm.gray)
#plt.title('blurred',size=20)
#plt.axis('off')
#plt.show()



#比较convolve和fftconvolve卷积函数的运行时间
#im_blurred2=signal.convolve(img,guass,mode='same')                            
#def wrapper_convole(func):
#    def wrapper_convole():
#        return func(img,guass,mode="same")
#    return wrapper_convole
#调用函数
#wra=wrapper_convole(signal.convolve)
#wra1=wrapper_convole(signal.fftconvolve)
#time1=timeit.repeat(wra,number=1,repeat=100)
#time2=timeit.repeat(wra1,number=1,repeat=100)
#plt.figure(figsize=(15,5))
#plt.gray()
#plt.subplot(131),plt.imshow(img),plt.title('Orig',size=15),plt.axis('off')
#plt.subplot(132),plt.imshow(im_blurred),plt.title('convel',size=15),plt.axis('off')
#plt.subplot(133),plt.imshow(im_blurred2),plt.title('fftconvel',size=15),plt.axis('off')

#绘制时间箱型图
#data=[time1,time2]
#plt.figure(figsize=(8,6))
#box=plt.boxplot(data,patch_artist=True)
#colors=['cyan','pink']
#for i,color1 in zip(box['boxes'],colors):
#    i.set_facecolor(color1)
#plt.xticks(np.arange(3),('','convole','fftconvolve'),size=15)
#plt.yticks(fontsize=15)
#plt.xlabel('scipy.signal meethods',size=15)
#plt.ylabel('time',size=15)
#plt.show()


#3.2频域滤波
#高通滤波
#读取图片并转为灰度图
#xn=np.array(Image.open('F:\\python练习\\图像实战练习\\images\\rhino.jpg').convert('L'))
##绘制原始图像的频谱
#freq=fp.fft2(xn)
#(w,h)=freq.shape
#half_w,half_h=int(w/2),int(w/2)
#freq1=np.copy(freq)
#freq2=fftpack.fftshift(freq1)
#plt.figure(figsize=(10,10)),plt.imshow((20*np.log10(0.1+freq2)).astype(int)),plt.show()
#
##应用HPF来阻断Numpy2D数组中的低频
#freq2[half_w-10:half_w+11,half_h-10:half_h+11]=0
#plt.figure(figsize=(10,10))
#plt.imshow((20*np.log10(0.1+freq2)).astype(int))
#plt.show()
#
##应用ifft2函数从上面的光谱中获取图像
#xn_img=np.clip(fftpack.ifft2(fftpack.ifftshift(freq2)).real,0,255)
#plt.imshow(xn_img,cmap='gray'),plt.axis('off'),plt.show()



#sy=np.array(Image.open('F:\\python练习\\图像实战练习\\images\\cameraman.jpg').convert('L'))
#快速傅里叶变换
#sy1=fp.fft2(sy)
#(w,h)=sy1.shape
#half_w,half_h=int(w/2),int(h/2)
#shrnp=[]
#lbs=list(range(1,25))
def signaltonoise(a, axis=0, ddof=0):
    a = np.asanyarray(a)
    m = a.mean(axis)
    sd = a.std(axis=axis, ddof=ddof)
    return np.where(sd == 0, 0, m/sd)
#plt.figure(figsize=(12,20))
##scipy已经移除这个函数
#for l in lbs:
#    sy2=np.copy(sy1)
#    sy3=fftpack.fftshift(sy2)
#    sy3[half_w-1:half_w+1+l,half_h-1:half_h+1+l]=0
#    sy4=np.clip(fp.ifft2(fftpack.ifftshift(sy3)).real,0,255)
##    print(signaltonoise(sy4,axis=None))
#    shrnp.append(signaltonoise(sy4,axis=None))
#    plt.subplot(6,4,l),plt.imshow(sy4,cmap='gray'),plt.axis('off')
#    plt.title('F='+str(1+l),size=20)
#plt.subplots_adjust(wspace=0.1,hspace=0)
#plt.show()


#绘制信噪比随截止频率变化的情况
#plt.plot(lbs,shrnp,'b.-')
#plt.xlabel('Cutoff',size=20)
#plt.ylabel('SNR',size=20)
#plt.show()



#低通滤波器
#fig,(axes1,axes2)=plt.subplots(1,2,figsize=(15,10))
#plt.gray()
#st=np.mean(io.imread('F:\\python练习\\图像实战练习\\images\\lena.jpg'),axis=2)
#st1=fp.fft2(st)
#st_gaussian=ndimage.fourier_gaussian(st1,sigma=4)
#st2=fp.ifft2(st_gaussian)
#axes1.imshow(st),axes1.axis('off'),axes2.imshow(st2.real)
#axes2.axis('off')
#plt.show()
#
##显示图像的频谱
#plt.figure(figsize=(8,8))
#plt.imshow((20*np.log10(0.1+fp.fftshift(st_gaussian))).astype(int))

#使用scipy中的fftpack实现低通滤波器
#bison=np.array(Image.open('F:\\python练习\\图像实战练习\\images\\bisons.png').convert('L'))
#bis_freq=fp.fft2(bison)
#(m,q)=bis_freq.shape
#half_m,half_q=int(m/2),int(q/2)
#bis_freq1=np.copy(bis_freq)
#bis_freq2=fftpack.fftshift(bis_freq1)
#bis_freq2_low=np.copy(bis_freq2)
#bis_freq2_low[half_m-10:half_m+11,half_q-10:half_q+11]=0
#bis_freq2-=bis_freq2_low
#bi_im=fp.ifft2(fftpack.ifftshift(bis_freq2)).real
#plt.imshow(bi_im,cmap='gray'),plt.axis('off')
#plt.show()


#绘制频谱（低频）
#plt.figure(figsize=(10,10))
#plt.imshow((20*np.log10(0.1+bis_freq2)).astype(int))
#plt.show()


#LPF作用在兔子身上的的应用
#bun=np.array(Image.open('F:\\python练习\\图像实战练习\\images\\bunny.png').convert('L'))
#bun_freq=fp.fft2(bun)
#(b,c)=bun_freq.shape
#half_b,half_c=int(b/2),int(c/2)
#bun_lp=[]
#ubs=list(range(1,25))
#plt.figure(figsize=(12,20))
#for u in ubs:
#    bun_freq1=np.copy(bun_freq)
#    bun_freq2=fftpack.fftshift(bun_freq1)
#    bun_freq_low=np.copy(bun_freq2)
#    bun_freq_low[half_b-u:half_b+u+1,half_c-u:half_c+u+1]=0
#    bun_freq2-=bun_freq_low
#    bun_im=fp.ifft2(fftpack.ifftshift(bun_freq2)).real
#    bun_lp.append(signaltonoise(bun_im,axis=None))
#    plt.subplot(6,4,u),plt.imshow(bun_im,cmap='gray'),plt.axis('off')
#    plt.title('F='+str(u),size=20)
#plt.subplots_adjust(wspace=0.1,hspace=0)
#plt.show()
  
 
#信噪比与截止频率变化关系
#bun_snr=signaltonoise(bun,axis=None)
#plt.plot(ubs,bun_lp,'b.-')
#plt.plot(range(25),[bun_snr]*25,'r-')
#plt.xlabel('Cutoff',size=20)
#plt.ylabel('SNR',size=20)
#plt.show()
    
   
#DOG带通滤波器
#tiger=img_as_float(io.imread('F:\\python练习\\图像实战练习\\images\\tiger.jpg'))
#plt.figure(),plt.imshow(tiger),plt.axis('off'),plt.show()
#x=np.linspace(-10,10,15)
#kernel_ld=np.exp(-0.005*x**2)
#kernel_ld /=np.trapz(kernel_ld)
#tiger_guass=kernel_ld[:,np.newaxis]*kernel_ld[np.newaxis,:]
#kernel_ld=np.exp(-5*x**2)
#kernel_ld /=np.trapz(kernel_ld)
#tiger_guass1=kernel_ld[:,np.newaxis]*kernel_ld[np.newaxis,:]
#DoGkernel=tiger_guass[:,:,np.newaxis]-tiger_guass1[:,:,np.newaxis]    
#tiger_im=signal.fftconvolve(tiger,DoGkernel,mode='same')
#plt.figure(),plt.imshow(np.clip(tiger_im,0,1)),print(np.max(tiger_im)),
#plt.show()    
##  


#3.2.5带阻滤波器
#plt.figure(figsize=(10,8))
#yu_im=np.mean(io.imread('F:\\python练习\\图像实战练习\\images\\parrot.png'),axis=2)/255
#plt.subplot(2,2,1),plt.imshow(yu_im,cmap='gray'),plt.axis('off')
#plt.title('Origal')
#F1=fftpack.fft2((yu_im).astype(float))
#F2=fftpack.fftshift(F1)
#plt.subplot(2,2,2),plt.imshow((20*np.log10(0.1+F2)).astype(int),
#cmap=plt.cm.gray)
#plt.xticks(np.arange(0,yu_im.shape[1],20))
#plt.yticks(np.arange(0,yu_im.shape[0],20))
#plt.title('Original')
#for n in range(yu_im.shape[1]):
#    yu_im[:,n]+=np.cos(0.1*np.pi*n)
#plt.subplot(2,2,3),plt.imshow(yu_im,cmap='gray'),plt.axis('off')
#plt.title('Image after adding noise')
#F1=fftpack.fft2((yu_im).astype(float))
#F2=fftpack.fftshift(F1)
#plt.subplot(2,2,4),plt.imshow((20*np.log10(0.1+F2)).astype(int),cmap=plt.cm.gray)
#plt.xticks(np.arange(0,yu_im.shape[1],20))
#plt.yticks(np.arange(0,yu_im.shape[0],20))
#plt.title('Noise')
#plt.tight_layout()
#plt.show()


#设计一个带阻或陷波滤波器
#F2[170:176,:220]=F2[170:176,230:]=0
#yu_im1=fftpack.ifft2(fftpack.ifftshift(F2)).real
#plt.axis('off'),plt.imshow(yu_im1,cmap='gray'),plt.show()


#3.2.6图像复原
#利用FFT去卷积和逆滤波
#len_im=255*color.rgb2gray(io.imread('F:\\python练习\\图像实战练习\\images\\lena.jpg'))
#len_gauss=np.outer(signal.gaussian(len_im.shape[0],3),signal.gaussian(len_im.shape[1],3))
#len_freq=fp.fft2(len_im)
#freq_kernel=fp.fft2(fp.ifftshift(len_gauss))
#convoled=len_freq*freq_kernel
#    
##模糊处理
#len_blur=fp.ifft2(convoled).real
#len_blur=255*len_blur/np.max(len_blur)
#
##对模糊图像使用逆滤波器
#epsilon=10**-6
#len_freq=fp.fft2(len_blur)
#freq_kernel=1/(epsilon+freq_kernel)
#convoled=len_freq*freq_kernel
#im_restored=fp.ifft2(convoled).real
#im_restored=255*im_restored/np.max(im_restored)
#print(np.max(len_im),np.max(im_restored))
#plt.figure(figsize=(10,10))
#plt.gray()
#plt.subplot(221),plt.imshow(len_im),plt.title('Original'),plt.axis('off')
#plt.subplot(222),plt.imshow(len_blur),plt.title('Blurred'),plt.axis('off')
#plt.subplot(223),plt.imshow(im_restored),plt.title('Restored'),plt.axis('off')
#plt.subplot(224),plt.imshow(im_restored-len_im),plt.title('Diffrence'),plt.axis('off')
#plt.show()


#利用维纳滤波器去卷积
#ele_im=color.rgb2gray(io.imread('F:\\python练习\\图像实战练习\\images\\whale.png'))
#from scipy.signal import convolve2d as conv2
#n=7
#psf=np.ones((n,n))/n**2
#ele_im1=conv2(ele_im,psf,'same')
#astro=img_as_float(ele_im1)
#ele_im1+=0.1* astro.std()*np.random.standard_normal(ele_im.shape)
#ele_im2,_=restoration.unsupervised_wiener(ele_im1,psf)
#fig,axes=plt.subplots(nrows=1,ncols=3,figsize=(10,4),sharex=True,sharey=True)
#plt.gray()
#axes[0].imshow(ele_im),axes[0].axis('off'),axes[0].set_title('Original',size=20)
#axes[1].imshow(ele_im1),axes[1].axis('off'),axes[1].set_title('Noisy',size=20)
#axes[2].imshow(ele_im2),axes[2].axis('off'),axes[2].set_title('set',size=20)   
#fig.tight_layout()
#plt.show()    
#    
    

#利用FFT实现图像去噪----灰度图像
#from matplotlib.colors import LogNorm
#moon_im=io.imread('F:\\python练习\\图像实战练习\\images\\moonlanding.png').astype(float)
#plt.figure(figsize=(10,10))
#plt.imshow(moon_im,plt.cm.gray),plt.axis('off'),plt.title('Orig'),plt.show()   
#
#
##傅里叶变换频谱
#moon_fft=fftpack.fft2(moon_im)
#def plot_spectrum(moon_fft):
#    plt.figure(figsize=(10,10))
#    plt.imshow(np.abs(moon_fft),norm=LogNorm(vmin=5),cmap=plt.cm.afmhot),plt.colorbar()
#plt.figure(),plot_spectrum(fftpack.fftshift(moon_fft))
#plt.title('Spectrum',size=20)
#
#
##FFT的滤波器
#keep_fraction=0.1
#moon_fft2=moon_fft.copy()
#r,c=moon_fft2.shape
#moon_fft2[int(r*keep_fraction):int(r*(1-keep_fraction))]=0
#moon_fft2[:,int(c*keep_fraction):int(c*(1-keep_fraction))]=0
#plt.figure(),plot_spectrum(fftpack.fftshift(moon_fft2)),plt.title('Filted')
##
#
##重建图像
#im_new=fp.ifft2(moon_fft2).real
#plt.figure(figsize=(10,10)),plt.imshow(im_new,plt.cm.gray)
#plt.axis('off')
#plt.title('Recon',size=20)
#


